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Abstract. We investigate the question of whether the recent modulation signal 
claimed by CoGeNT is best explained by the dark mattcir (DM) hypothesis from 
a Bayesian model comparison perspective. We consider five phenomenological 
explanations for the data: no modulation signal, modulation due to DM, modulation 
due to DM compatible with the total CoGeNT rate, and a signal coming from other 
physics with a free phase but annual period, or with a free phase and a free period. In 
each scenario, we assign to the free parameters physically motivated priors. We find 
that while the DM models are weakly preferred to the no modulation model, but when 
compared to models where the modulation is due to other physics, the DM hypothesis 
is favoured with odds ranging from 185 : 1 to 560 : 1. This result is robust even when 
astrophysical uncertainties are taken into account and the impact of priors assessed. 
Interestingly, the odds for the DM model in which the modulation signal is compatible 
with the total rate against a DM model in which this prior is not implemented is only 
5 : 8, in spite of the former's prediction of a modulation amplitude in the energy range 
0.9 — )• 3.0 keVee that is significantly smaller than the value observed by CoGeNT. 
Classical hypothesis testing also rules out the null hypothesis of no modulation at the 
1.6(7 to 2.3(7 level, depending on the details of the alternative. Lastly, we investigate 
whether anisotropic velocity distributions can help to mitigate the tension between the 
CoGeNT total and modulated rates, and find encouraging results. 
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1. Introduction 

As suggested by a host of independent cosmological and astrophysical observations, 
most of the matter content in the Universe is invisible and reveals its existence only 
through gravitational effects. This is the so-called dark matter (DM). Amongst the most 
well-motivated candidates to explain this dark matter are Weakly Interacting Massive 
Particles (WIMPs), which have weak-scale cross-sections and masses in the GeV to 
TeV range. Many varied efforts are currently underway to search for these particles. 
In particular, direct detection experiments have been designed to observe WIMPs that 
have clustered in the Galactic halo, through their coherent scattering with nuclei in a 
target material. 

The expected WIMP scattering rates in direct searches are generally small, and the 
exponential decay of the WIMP recoil-spectrum mimics that of the background at low 
energies. However, a specific signature can be exploited in order to disentangle a DM 
signal from the background, namely, the annual modulation of the scattering rate [1, 2]. 
This effect, expected to be of the order of a few percent, is induced by the rotation 
of the Earth around the Sun, which periodically changes the DM flux relative to the 
detector, and is known to be a smoking gun signature for DM detection. 

Observation of such an annual modulation signal has long been claimed by the 
DAMA/LIBRA experiment, a detection that has by now reached a significance of 
8.9a [3,4], and can be interpreted in terms of the scattering of a light WIMP of mass 
~ 10 GeV off sodium nuclei. Interestingly, a similar light WIMP explanation appears 
to explain also the low energy excess reported by the CoGeNT experiment in 2010 [5] , 
as well as the more recent CRESST results [6]. Intriguingly however, these results are 
at the same time excluded by a number of null searches, most notably XENON [7], and 
upon close scrutiny, are also mildly inconsistent with one another [8-10]. Adding further 
to the puzzle, the CoGeNT collaboration recently reported a modulating signal in the 
energy range 0.5 — ?> 3.0 keVee at 2.8a, with a modulation amplitude of 16.6 ± 3.8%, 
a period of 347 ± 29 days, and the minimum rate falling on October 16±12 days [11], 
somewhat at odds with the DM prediction of December 2. 

The CoGeNT modulation result has sparked a new frenzy of searches for particle 
physics solutions that could reconcile all existing experimental findings, e.g., [12-30], as 
well as astrophysical solutions based on, e.g., streams [31] and anisotropic DM velocity 
distributions [32]. In this work, we examine the CoGeNT modulation data from a 
Bayesian model comparison perspective [33-35]. We ask the question of whether the 
data really show the presence of a modulation, and if so, whether this modulation is 
best explained in terms of WIMP scattering or by some other physical process. 

The classical approach of hypothesis testing does not and cannot give the 
probability for a hypothesis, which is not defined in frequentist statistics. Rather, it 
returns the probability of observing as extreme or more extreme values of the test 
statistic assuming the null hypothesis is true, i.e., the p- value. This, however, is in 
general not the scientific question one is interested in (see, e.g., [35, 36] for a review). In 
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order to explain a phenomenon, one would like to assess the probability of competing 
models under the data. This can be done only in the context of Bayesian inference, where 
a probability can be assigned to all propositions — be they propositions of a parameter 
value within a given model (such as in the practice of Bayesian parameter inference), 
or of the model itself. In this way, the task of selecting the "best" amongst competing 
models to describe a set of data becomes formally one of finding the model whose 
posterior probability is the highest. 

The central quantity in Bayesian model comparison is the evidence, or marginal 
hkehhood, of the model, i.e., the hkelihood function for the model's parameters 
averaged over parameter space weighted by the prior probability of the parameters. 
By construction the evidence automatically incorporates the notion of Occam's razor; 
models that fit the data well are rewarded through a favourable likelihood function; 
models that are unpredictive (e.g., excessively broad priors) are penalised by the 
larger parameter volume over which the likelihood must be averaged. Bayesian model 
comparison has found widespread application in cosmological data analysis, from 
curvature testing [37] to infiationary model selection [38,39], amongst others [40-42]. 
Its use in particle physics is less common, but see, e.g., [43-45]. 

Lastly, although the main thrust of our analysis is Bayesian in nature, we also 
highlight along the way some of the common pitfalls in the computation of classical 
p-values, especially in the context of DM direct searches. We use this opportunity to 
recalculate the p-values for DM hypotheses that have been misestimated in previous 
works. 

The rest of the article is organised as follows. We introduce the statistical 
underpinnings of our analysis in section 2, and describe the analysis of the CoGeNT 
data in section 3. After a brief review of the theory of modulation due to WIMP 
scattering in section 4, we define the models for testing in section 5. Our results are 
presented in section 6. Section 7 contains our conclusions. 

2. Statistical approach 

2.1. Bayesian inference 

Given a set of parameters 6 defining a model Ai, we are interested to compute their 
posterior probability distribution function (pdf) p{9\d,A4) via Bayes' theorem, namely. 



Here, d arc the data under consideration, C{6\Ai) = p{d\6,Ai) the likelihood function, 
and p{9\Ai) is the prior pdf for the parameters under the model. The quantity p{d\M), 
defined as 



p{9\d,M) 



c{e\M)p{e\M) 

pid\M) 





(2.2) 



is called the Bayesian evidence. 
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Bayesian inference is based on the posterior pdf for the parameters 6, and it assumes 
that the model under consideration, A4, is the correct one. But we can expand our 
inferential framework to ask the question of the viability of the model itself, or rather, 
of the relative performance of alternative possible models as explanations for the data 
— this is the subject of Bayesian model comparison. The formalism of Bayesian model 
comparison automatically balances the quahty of the model's fit to the data against 
its predictiveness. Thus, if there are several competing models, the problem of finding 
the "best" model — one that achieves the optimum compromise between quality of fit 
and predictiveness — can be formally defined as selecting the model that has the highest 
posterior probability. In this sense, the methodology of Bayesian model selection can 
be interpreted as a quantitative expression of the Occam's razor principle of simplicity. 

We now turn to a brief review of the framework of Bayesian model comparison. For 
a more in-depth discussion see, e.g. [34,35]. 



2.2. B ayes factors and model comparison 

From the definition of the Bayesian evidence in equation (2.2), it is easy to see that 
this quantity incorporates the notion of Occam's razor and penalises those models with 
excessive complexity unsupported by the data for wasted parameter space. Increasing 
the dimensionality of the parameter space without significantly enhancing the likelihood 
jC{d\9,A4) in the new parameter directions reduces the evidence. Unpredictive priors 
p{9\M.) (e.g., excessively broad compared with the width of the likelihood) hkewise 
dilute the evidence. 

The posterior probability p{Ai\d) of a model A4 given the data d is related to the 
Bayesian evidence via Bayes' theorem, 

p{M\d) (X p{d\M) p{M) , (2.3) 

where p{M) is the prior probabihty assigned to the model M. itself, and we have dropped 
a normalisation constant corresponding to the probability of the data d. Thus, the 
posterior odds between two competing models A4o and A4i are given by 

p{M,\d) R P(-Mi) 

~ -Pio / . > X , (2.4) 



where 



p{Mo\d) 'piMo) ' 



_ p{d\Mi) 

= pidlMo) ' ^'-'^ 
a ratio of the models' evidences, is called the Bayes factor. 

The Bayes factor Biq represents an update from our prior belief in the odds of two 
competing models p{/Ai) / p{/Ao) to the posterior odds p{A4i\d) / p{Aio\d) ■ If there are 
no prior reasons to believe that one model should be more probable than the other, 
then p{Aii) = p(A^o) (non-committal prior), and the Bayes factor alone determines the 
outcome of the model comparison. A Bayes factor larger than unity means that the 
model All is preferred over the model A^o ^ a description of the experimental data. 
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Table 1. Jeffreys' scale for grading the strength of evidence for two competing models 
M.0 and A4i, adapted from [35,38]. 



In Bio 


Odds Ml : Mo 


Strength of evidence 


< -5.0 


< 1 : 150 


Strong evidence for Mo 


-5.0 ^ -2.5 


1 : 150 ^ 1 : 12 


Moderate evidence for Mo 


-2.5 ^ -1.0 


1 : 12 ^ 1 : 3 


Weak evidence for Mo 


-1.0 ^ 1.0 


1:3-^3:1 


Inconclusive 


1.0 ^ 2.5 


3 : 1 ^ 12 : 1 


Weak evidence against Mo 


2.5 5.0 


12 : 1 150 : 1 


Moderate evidence against Mo 


> 5.0 


> 150 : 1 


Strong evidence against Mo 



and vice versa. As may be expected, the correspondence between the actual value of 
the Bayes factor and the degree of preference in the ordinary sense — and hence the use 
of the Bayes factor as a decision-making criterion — is a matter of convention.! Here, we 
use the convention set down by "Jeffreys' scale" shown in table 1. 

As can be seen in equation (2.2), the computation of the evidence p{d\M) for each 
model M requires the evaluation of an integral over the parameter space. For this 
purpose we use the multimodal nested-sampling algorithm implemented in the publicly 
available package MultiNest v2.12 [46,47]. We set niive = 10000 and a tolerance factor 
of 0.01, following the recommendations of [46]. 

2.3. Sensitivity analysis for the Bayes factors 

The evidence and the Bayes factors are obviously dependent on the choice of the prior 
probability distributions for the model's parameters, p{6\M). Since the choice of priors 
is usually not unique, interpretation of the results of Bayesian model selection ought 
to allow for the impact of a reasonable change of priors. This is called "sensitivity 
analysis" . 

As we shall describe in detail in section 5, our approach to studying time modulation 
in the CoGeNT data is purely phenomenological. We begin with the parameter space 
of the most complex model (in terms of the number of free parameters it contains), 
^ ~ {^m' 'S'm) ^max? ^}) ^nd define from it a set of nested models, where simpler models 
are realised by setting some or all of these parameters to specific values 9*. With one 
exception (model lb, sec table 5), the models we test are not tightly connected to the 
predictions of any specific physical process. This approach, while perfectly convenient 
for the purpose of parameter inference, poses a problem for model comparison: without 
the guidance of specific predictions, the odds for a more complex model can be made 

I In the same way, the classical statistics criterion that a null hypothesis should be rejected if the 
p-value falls below 0.05 is a matter of convention. 
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arbitrarily small by increasing the width of the priors on the additional parameters or 
by choosing uniform priors on non-linear functions of these parameters. 

If the models Aio and Aii are nested and their parameter priors separable, then 
the impact of changing the prior width on the Bayes factor can be estimated analytically 
using the Savage-Dickey density ratio (SDDR, see [33]), 

Introducing the notation 9 = {^,ip), where denotes the additional parameters of 
the more complex model Aii, and ip labels the free parameters common to both models, 
p{i!}*\d, Ail) = A1i)d'0 is the marginal posterior pdf for the additional 

parameters of A4i evaluated at A^o's parameter value marginalised over ip, and 
^(■j^^lA^i) = ■0|A1i)d'0 is the marginal prior density of Ali evaluated at the 

same 

If the prior p(^9|A^i) is sufficiently broad and the data sufficiently constraining, 
then p{'d*\d, Ail) will be approximately independent of the prior. However, because 
the prior pdf must be normalised to unity probability content, increasing the width 
of p('i?|A^i) leads to a smaller p{'i)*\Aii) and therefore to a smaller Bayes factor. If 
for instance p{-d\Aii) is a top-hat function, the SDDR formula shows that rescahng its 
width by a factor a will change In Biq by approximately — In a. We will use this analytic 
approximation to perform a sensitivity analysis of our model comparison results. 



3. CoGeNT data analysis 

The CoGeNT collaboration has recently released data corresponding to 458 days of data 
taking, of which 442 days are live [11]. The initial day tin is December 4, 2009, from 
which point on data taking has been continuous until March 6, 2011, except for the 
gaps on days 68 — > 74, 102 — > 107, and 306 — > 308 inclusive. These gaps are taken into 
account in our time-binning of the data. The fiducial detector mass is 330 g, leading 
to a total exposure of 145.86 kg-days. Although the energy threshold of the detector is 
nominally 0.4 kcVcc (electron equivalent keV), we consider in this work a threshold of 
0.5 keVee in order to avoid trigger threshold effects. 

The signal region receives contributions from events induced by cosmogenic 
activation of radioisotope decays via electron capture. Several peaks arc expected 
at energies predicted by i^T-shell decay ranging from 0.56 to 1.4 keVee, the dominant 
being the Ge^^ and the Zn^^ peaks at 1.1 and 1.3 keVee, respectively. These can be 
subtracted following the prescription of the GoGeNT collaboration [48]. Henceforth, we 
shall present event rates only after subtracting these contributions. 

Our main focus in this work is the time modulation of the CoGeNT data. It 
is on the measured event rate as a function of time that we perform our Bayesian 
model comparison analysis. However, as we shall see in section 5, the total number 
of events summed over the whole measurement period can be used to construct self- 



Evidence for dark matter modulation in CoGeNT? 



7 




100 200 300 400 500 100 200 300 400 500 100 200 300 400 500 

t (days) t (days) t (days) 



Figure 1. CoGeNT time-stamped data binned in intervals of 30 days, in three energy 
bins. From left to right, AE^ ^ 0.5 0.9 keVee, AE2 = 0.9 3.0 keVee, and 
AE^ = 3.0 —7- 4.5 keVee. Note that contributions from radiative peaks have been 
removed. In each panel, the black curve denotes the best-fit point for model la, the 
red dot-dashed line for lb, the green long-dashed line for 2a, while the dotted gray line 
indicates the constant component of the event rate. 



consistent priors on the modulation amplitude. We describe the analysis of both the 
time- dependent event rate and the total event rate below. 

3.1. Time- dependent rate 

We investigate the time evolution of the event rate in three top-hat energy bins (index i): 
AEi = 0.5 ^ 0.9 keVee, AE2 = 0.9 ^ 3.0 keVee, and AE3 = 3.0 ^ 4.5 keVee. 
The choice of these bins is based on theoretical expectations for a light-mass WIMP 
interacting in a CoGeNT-like detector: there should be a time modulation in the signals 
in z = 1, 2, and none in i=3. For each energy bin, we group the time-stamped data in 
intervals of 30 days starting from tin. The last interval of 8 days is discarded from the 
analysis, because of its minor statistical weight due to the large error bar. This gives a 
total of 15 time bins (index j). We show in figure 1 the time-binned data Cij for the 
three energy bins, in units of counts per 30 days. The corresponding error bars (jjj are 
estimated from Poisson statistics, i.e., CTjj = {30 /xj)\l Cij, where Cij is the raw number 
of events prior to subtraction of the radiative peaks in the ith energy bin and the jth 
30-day interval, and xj is the actual number of days of data taking in the jth time 
interval, in the event of gaps. 

The likelihood function for the ith energy bin is approximated as an uncorrelated 
multivariate Gaussian function, i.e., up to irrelevant normalisation constants, 

j=l ^^ij 

where Rij denotes the theoretical event rate in the ith energy bin and jth time bin. 
In the case of WIMP scattering, Rij is a function of the DM parameters as well as 
astrophysical quantities. We discuss in more detail the computation of Ri^ in section 5. 
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3.2. Total rate 

We do not use the CoGeNT total event rate directly in our model comparison analysis. 
However, as we shall explain in section 5, in some cases we are interested to estimate 
the DM parameters preferred by the total rate, and use this information to construct 
self-consistent priors on the modulation amplitude to be used in the model comparison 
analysis of the time-dependent data. 

To this end, we group the total number of events observed in the energy range 
0.5 — )■ 3.2 keVee into 27 energy bins of width = 0.1 keVee. In each of these bins, 
the quantity Ck denotes the number of counts per day per detector mass per energy (in 
units cpd/kg/keVee) over the period of data taking, with the corresponding statistical 
error ak estimated from Poisson statistics prior to subtraction of the radiative peaks. 
The likelihood function for the total rate is approximated as a Gaussian, i.e., 

InAotai^-E ^^^^+M"^^^' , (3.2) 
k=i ^^fe 

where is the predicted DM signal rate in the kth energy bin, while is the background 
contribution to that bin. Note that by approximating the likelihood function as an 
uncorrelated multivariate Gaussian, we have implicitly assumed the count rates in the 
individual energy bins to be uncorrelated. This should be a reasonable assumption. 



given that the CoGeNT energy resolution, ae ~ + 981£^/keVee keVee [49], i.e., 

a£ ~ 0.02 — )■ 0.057 keVee in the energy range considered in this work, is always smaller 
than the energy bin width A5. 

The DM signal Sk{mY)M, o'^', Pq, f{v')) comes from integrating the differential recoil 
rate for coherent and elastic WIMP scattering off nuclei over the energy range of the 
kth bin, and is a function of the DM mass moM and the spin-independent cross-section 
(T^^, as well as of the local DM density Pq and several other astrophysical parameters 
characterising the DM velocity distribution f{v') in the solar neighbourhood. We shall 
elaborate more on the astrophysics in section 5 where appropriate. For details of the 
computation of Sk, we refer the reader to [50]. Note however a technical detail that 
differs in the present analysis: the form of the (energy-dependent) quenching factor 
qcc, which relates the observed ionisation energy S (in units keVee) to the actual recoil 
energy E (in units of nuclear recoil keV, keVnr) via S = qceE. Here, we use Lindhard's 
theory with k — 0.2, i.e., 

^(keVee) = 0.19935 x ^^•^^''^(keVnr) , (3.3) 

as recommended in [48]. 

The background rate is obtained similarly from integrating over the A:th bin the 
differential background dB/d£, which we model as a constant factor plus an exponential 
decay, 

^=C + ^exp(-£/£o), (3.4) 



Evidence for dark matter modulation in CoGeNT? 



9 



Table 2. Priors for the nuisance parameters describing the background in the CoGeNT 
experiment. All priors are uniform over the indicated range. 



Parameter 


Prior 


So 


^ 30 (kcVcc) 


C 


^ 10 (cpd/kg/keVee) 


A 


^ 10 (cpd/kg/keVee) 



so that 



1 

C + 



dS 



exp - 



— exp 



gfc + Afc 
So 



(3.5) 



with nuisance parameters C, A and Eq. The exponential background has been introduced 
to account for surface event contamination in the signal region due to misidentified 
electrons in the nuclear recoil band. This and the flat background rate C are assumed 
to be constant in time. Table 2 shows the prior ranges for these nuisance parameters. 



4. Annual modulation due to dark matter 



Direct dark matter searches are sensitive to both the particle physics of the DM 

— * 

candidate and the local properties of the Galactic halo in the solar neighbourhood, Rq. 
The event rate due to scattering of WIMPs in a dark matter detector is proportional to 

ri{E,t)=l dh'i^, (4.1) 

where f{v'{t)) is the local velocity distribution of the DM particles with respect to the 
Earth's frame (primed), and Vmin is the minimum velocity required of the particles to 
deposit an energy E in the detector. See for further details [50] and references therein. 

To compute f{v'{t)), we note that on the timescale of a direct detection experiment, 
the local DM velocity distribution in the Galactic frame (unprimed), F{v,Rq), can be 
taken to be constant in time. However, because the solar system moves through the DM 
halo with a velocity Vq in the Galactic frame, and the Earth moves around the Sun with 
a rotation velocity f "©,rot in the Sun's rest frame, the local DM velocity distribution as 
seen on Earth, f{v'{t)), becomes time-dependent via the transformations 

^^e = |^^0 + ^^"e,rot| = + ^^e,rot cos 7 cos[27r(t -to)/T], (4.2) 

where 7 = 60° is the inclination of the Earth's rotation plane with respect to the Galactic 
plane, to ~ June 2, and T = 1 year. Note that the second equation is a projection of 
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the Earth's orbital motion onto the Galactic plane and neglects the small eUipticity of 
the orbit. 

We can easily estimate the effect of the Earth's movement on the DM scattering 
rate in the |f' + -U0| <^ v'^^.^^. limit by expanding ri(E,t) around its "mean" value r]Q{E,t) 
(i.e., evaluated at ^^rot — 0) terms of the small parameter e = v'^^^^^/\v' + Vq\. For 
an isotropic halo, i.e., F{y, R) = F{v, R), we find 

7]{E, t) - TjoiE, t) oc cos[27r(t - to)/T] , (4.3) 

which is the well-known result that the WIMP scattering rate has a sinusoidal time 
dependence with a period of one year, and a phase determined by the movement of the 
Earth with respect to the Galactic frame, %,. We shall deal only with isotropic DM 
velocity distributions in our model comparison analysis, i.e., equation (4.3) is always 
taken to be the prediction of WIMP scattering. We note however that anisotropic 
distributions, streams, or unvirialised components in the halo phase space may induce 
changes in the phase of modulation as well as in its functional form [31, 51-55]. 

5. Model definitions 

The main goal of this paper is to evaluate in a quantitative manner the probability of 
various explanations (i.e., models) for the presence (or absence) of modulation in the 
GoGeNT time-dependent data, and in particular to estimate the probabihty of the data 
being due to a hght-mass WIMP signal. To this end, we adopt a phenomenological 
approach, and, inspired by equation (4.3), parameterise the time-dependent event rate 
R{t) in the ith energy bin as 

Ri(t) = C/i (l + Si, cos[27r(i - W - 28)/r]) , (5.1) 

where t is in units of days, with t — corresponding to the first day of data-taking, i.e., 
December 4, 2009, t^ax is the phase of the modulation in terms of days since January 1, 
and T the modulation period. Two more parameters, and ^S"^, denote the mean event 
rate and the fractional modulation respectively, with the superscript "f indicating that 
these quantities are generally energy-dependent. The binned rates Rij are then simply 
given by 

Rij - ^ Riit) dt, (5.2) 

where At = 30 days. 

As discussed in section 4, if the scattering signal is due to dark matter in an isotropic 
halo, then the modulation phase and period are fixed at imax = 152 days (i.e., June 2) 
and T — 365 days respectively, while and S^^ are determined by the DM mass 
and cross-section as well as the background signal. In our phenomenological analysis, 
however, we treat all or subsets of {U^, Sj^, tmax- T} as free parameters, and determine 
using Bayesian model comparison if the CoGeNT data bear out the dark matter 
hypothesis or some other alternative scenarios. We describe the various hypotheses 
we wish to test below; a summary can be found at the end of the section in table 5. 
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Table 3. Priors for the mean rate in the three energy bins. All priors are uniform 
in the indicated ranges. 



Bin (keVee) Prior 

A^i = 0.5 ^ 0.9 20 60 (counts/30 days) 

AE2 = 0.9 ^ 3.0 30 ^ 100 (counts/30 days) 

AE3 = 3.0 4.5 10 60 (counts/30 days) 



5.1. Model 0: no modulation 

This is our baseline and also simplest model, representing the case of no modulation 
in the data, i.e., -S"^ = -S"^ = -S"^ = 0. In other words, the only relevant and free 
parameters are U^, and U^. These are allowed to vary in the ranges indicated in 
table 3, under uniform priors. The prior ranges have been chosen in such a way that they 
more than encompass the highest and the lowest measured count rates — including their 
error bars — in each energy bin. Note that although the prior ranges for these mean rate 
parameters are somewhat arbitrary, they do not impact on the outcome of our model 
comparison, because all models considered here have the same parameters and their 
prior volume cancels out by virtue of the Savage-Dickey density ratio formula (2.6) 
(see [33] for details). 

5.2. Model 1: modulation due to dark matter 

This set of models corresponds to modulation due to a DM signal, where the modulation 
period is one year, and the phase June 2 is predicted by the scattering of DM in an 
isotropic halo. We consider two specific implementations. 

Model la: phenomenological DM. Here we impose a uniform prior in the range — ?> 0.2 
on the fractional modulation amplitudes S}^ and S'j^ in the first two energy bins AEi and 
AE2, while for AE^ we assume no modulation, i.e., = 0. Such a model encodes the 
"naive" prediction for the modulated signal due to WIMP scattering in a CoGeNT-like 
detector [1, 2, 56]. Compared with model 0, model la contains one extra free parameter 
in each of AEi and AE2, while there is no change in AE3. 

Model lb: consistent DM signal. The choice of priors on Sl^ reflects our prior belief 
in the degree of modulation we expect to see in the experimental results under the 
DM hypothesis. This belief and hence the priors may be updated as we acquire more 
information from other experiments. In this sense, if our prior belief is that the CoGeNT 
experiment has indeed detected WIMP scattering events and that the modelling of the 
background is well understood, we can then use the CoGeNT total event rate to update 
the prior on the fractional modulation amplitude in each energy bin. 
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Table 4. Astrophysical constraints on the DM halo profile and the WIMP velocity 
distribution for model lb. Except for the halo concentration Cvir which is subject to 
a theoretical prior that is uniform in the indicated range, all other constraints are 
observational and in the form of Gaussians. 



Observa.l)le 



Constraint 



Local standard of rest 
Escape velocity 
Local DM density 
Virial mass 



vo = 230 ± 24.4 km s"^ [60, 61] 
^;esc = 544 ± 39 km s~^ [62, 63] 
Pq = 0.4 ± 0.2 GeV cm"^ [64, 65] 
Mvir = 2.7 ± 0.3 X 10^2 [66, 67] 



Halo concentration 



Cvir = 5 ^ 20 [68] 



To achieve this purpose, we use Markov Chain Monte Carlo (MCMC) techniques to 
infer the joint posterior pdf of the DM mass and cross-section — as well as the nuisance 
and astrophysical parameters — preferred by the CoGeNT total event rate. The DM 
mass is sampled using a uniform prior in the range 1 ^ 20 GeV, while a uniform prior 
on log(o"^Vcm^) in the range —41 — )■ —39 is assumed for the cross-section. For the DM 
velocity distribution, we consider an isotropic distribution arising from a Navarro-Frenk- 
White (NFW) density profile [57].§ Technically, the parametric NFW profile, with its 
free two parameters M^ir, the virial mass, and Cvk, the halo concentration, is "inverted" 
using the Eddington formula under the assumption of hydrostatic equilibrium. The 
resulting DM velocity distribution is subject to a number of observational and theoretical 
constraints listed in table 4. For more details see [50]. 

From the joint posterior pdf of {mDM, • • •} one can easily derive posterior pdfs 
for the fractional modulation amplitudes 5"^ by evaluating for every sample in the 
Markov Chain the expression 



energy bin, and June 2 (December 2) is the day when the WIMP scattering rate is 
expected to reach a maximum (minimum). Importantly, the background event rate, as 
given in equation (3.5), is treated as a constant component in time. The resulting chains 
in 5"^ are binned to construct the posterior pdfs, which are in turn the priors we seek 
for the present modulation analysis. Note that because we are using only the total rate 
to derive a prior for 5^ and no time information at this stage, we are not fitting the 
same data twice in the procedure. 

§ Although it is known that the Milky Way and spiral galaxies in general are better described by cored 
profiles [58,59], the choice of halo density profile has little impact on the inference results from the 
current generation of direct DM searches, as shown in [50]. 
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Figure 2. Priors for the fractional modulation amplitudes S"^ in model lb, derived 
from the CoGeNT total rate. 
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Figure 3. Samples from the pos- 
terior pdf in the {todm, cr^^, S*^.,} 
plane, where the S^^ value is indi- 
cated by the colour coding. The 
top panels show (left) and 5^ 
(right), while the bottom panel 
shows S*,?,. 



In the case we consider, inclusion of astrophysical nuisance parameters in deriving 
effective priors for leads to the width of such priors being larger than it would 
otherwise have been, had the astrophysical quantities been kept fixed. This in turn 
means that the evidence value for model lb will be reduced by virtue of the Occam's 
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razor effect with respect to the case in which the astrophysical parameters are held fixed. 
In this sense, our evidence values for model lb arc to be considered conservative. 

Figure 2 shows the priors on S'^j thus constructed for the three energy bins, while 
figure 3 shows the correlations between S*^, the DM mass, and the cross-section. In the 
case of the first bin, the fractional modulation amplitude increases with the cross- 
section, while its dependence on the DM mass is weak. In contrast, S"^ in the second 
bin depends strongly on the DM mass, with higher masses corresponding to larger 
modulation amplitudes. These observations indicate that and are minimally 
correlated with one another, and thereby justify our treatment of them as independent 
parameters in the model comparison analysis. Finally, in the third bin is close to 
zero, since for WIMPs with masses smaller than ~ 10 GeV, the exponentially decaying 
DM scattering signal becomes vanishingly small. 

5. 3. Model 2: modulation due to some other physics 

If the time modulation is not due to WIMP scattering off nuclei but to some unknown 
physics, then the phase tmax need not correspond to June 2, or the period T exactly one 
year. In this case, then, tmax and T should also be treated as free parameters alongside 
Sl^ and U^. We consider two distinct scenarios. 

Model 2a: non-DM, annual modulation. This model has an annual modulation, but 
with its modulation amplitudes and phase treated as free parameters. For the latter, we 
vary imax within — > 365 days under a uniform prior. Note that unlike the parameters 
Sl^ and Ul^ which can change from energy bin to energy bin, there is only one tmax 
parameter in the analysis. In other words, in a combined analysis of all three energy 
bins, the time modulations in the expected signals are described by the same phase. 

Model 2b: non-DM, free period. This is similar to model 2a, but in addition, the period 
T is allowed to vary in the range 1 365 days under a uniform prior. As with the 
parameter tmax, we use only one T parameter in the analysis. 

All models under consideration are summarised in table 5. 
6. Results and discussions 

The main results of our analysis are the Bayes factors, relative to model 0, for 
models 1 and 2. These are displayed in figure 4, and the corresponding odds for each 
model against model are listed in table 6. For reference, we show also the difference 
in twice the best-fit log-likelihood value, Ax^g, and, where analytically calculable, the 
p- values of the null in table 7. We elaborate on the sahent features of our results below. 
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Table 5. Summary of all models studied here. Unless otherwise stated, all priors 
are uniform in the indicated ranges. The quantity u counts the number of extra free 
parameters in a model with respect to model 0, with the first number referring to the 
one-bin analysis and the second to the all-bins analysis. 



Model 


Description 


Fractional 


Phase tynax 


Period T 


Extra 






modulation 


(days) 


(days) 


params 





No modulation 









z/ = 0,0 


la 


Pheno-DM 


^ ^ 0.2 

sl = o 


152 


365 


u = l,2 


lb 


Consistent DM 


Gaussian, clipped at 
(5^ > 0) 

Si = 0.098 ± 0.021 
= 0.026 ±0.011 
5^ = (0.37 ±36) X 10-^ 


152 


365 


z/ = 1,3 


2a 


Non-DM, annual 


0^1 


0^ 365 


365 


z/ = 2,4 


2b 


Non-DM, free period 


0^1 


0^ 365 


1 ^ 365 


v = 3,5 



6.1. Bin-by-bin description 

Energy bin 1 As shown in figure 4, the evidence for modulation in AE'i is generally 
inconclusive for DM models (models 1), and weakly against "other physics" models 
(models 2), when compared with the no modulation model. This yields a moderate 
to strong evidence for DM models as compared to the "other physics" models. Of the 
DM models, model lb (DM signal consistent with CoGeNT total rate) is marginally 
favoured over model la (phcnomenological DM model), despite the fact that purely 
from the quality of fit (see table 7), model lb performs marginally worse than la. This 
is an example of Lindley's paradox (i.e., Bayesian model selection returning a different 
result from classical hypothesis testing, see [33] and references therein): model lb is 
rewarded with a larger Bayes factor because of its narrow prior, S]^ — 0.098 ± 0.021, 
which makes it highly predictive compared with model la (prior S]^ = — )■ 0.20). 

The posterior marginal distribution for the fractional modulation amplitude in each 
scenario is shown in the top left panel of figure 5. Observe that all DM scenarios select 
comparable values for the modulation amplitude, e.g., model la prefers = 0.08±0.03 
and model lb S^^ = 0.09 ±0.01, because in both cases the phase has been fixed to June 2 
(^max ~ 152 days). However, if the phase is let free to vary, model 2a selects a slightly 
larger value for the modulation amplitude, 5"^ — 0.11 ± 0.04, and a preferred phase 
^max = 106 ± 29 days, as shown in figure 6. Model 2b is heavily multimodal both in t^ax 
and the period T (see also figure 7, top left panel). Because of this multimodahty, we 
do not summarise the inference for tmax ^-nd T in model 2b in terms of the ID posterior 
mean and standard deviation. 
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Figure 4. Bayes factors for the various modulation scenarios analysed in this work. 
The models are specified on the vertical axis, while the different symbols refer to the 
energy bin(s) for which the Bayes factors have been computed, as labelled in the plot. 
The actual value of InB in each case is indicated by the number above the data point. 
The Bayes factors have uncertainties of ~ 0.02 for the individual bins and ~ 0.04 for 
the combined analysis. Following Jeffrey's scale in table 1, the vertical lines demarcate 
the different empirical gradings of the strength of the evidence. 



Energy bin 2 Compared with model 0, the support for modulation in energy bin 2 
is at best weak. With the exception of the moderate evidence for model la against 
model 2b, the comparisons between the modulation models themselves are likewise weak 
to inconclusive. The preferred fractional modulation amplitude is = 0.20 ± 0.05 for 
both models 2a and 2b, while for model la the posterior mode is close to the prior 
boundary of 0.2 (see figure 5). By comparing the posterior for S"^ for models 2a and 
2b with that for models la and lb, it is clear that the data in energy bin 2 prefer a 
larger modulation amplitude that can be provided by either DM models. The posterior 
pdfs for both models 2a and 2b are unimodal in tma.x and, for model 2b, T, giving 
i^max = 104 ± 11 days for model 2a, and tmax = 108 ± 15 days and T = 331 ± 23 days for 
model 2b, as shown in figures 6 and 7. 
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Figure 5. ID marginal posterior pdfs (normalised to the peak) for, from left to right, 
the modulation amplitudes S^, and in the 3 CoGeNT energy bins, for all 
models considered in this work. Results from the single bin analyses are displayed 
in the top panels, while the bottom panels show results from the combined analysis, 
which includes all 3 energy bins in the likelihood. 
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Figure 6. ID marginal posterior pdfs for the phase imax in models 2a and 2b (left 
and middle panel respectively), and and for the period T in model 2b (right panel) 
from the single bin and the combined analyses, as labelled in the plots. 



Energy bin 3 For WIMP masses smaller than 10 GeV, the generic prediction for 
CoGeNT-like detectors is that there should be little to no annual modulation in the 
3.0 — )■ 4.5 keVee energy range. Indeed, we find in this energy bin that the no modulation 
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scenario is au pair with the DM models, and weakly to moderately preferred over the 
modulation due to "other physics" models. It is interesting to note that the evidence 
against model 2b relative to model is only weak. This is because with its freely varying 
phase tmax and period T, model 2b turns out to provide a better fit to the time-structure 
of the CoGeNT data (left panel of figure 1) than a constant rate. As can be seen in 
figures 6 and 7, the posterior pdf is unimodal in the period T, with a preference for 
T = 91 ± 4 days, but is strongly multimodal in the direction of the phase tmax- 

Thus, we conclude that in the energy range 3.0 — > 4.5 keVee, the CoGeNT data do 
not support the presence of modulation. This statement is reinforced by the behaviour 
of the ID marginal posterior for in the top right panel of figure 5, which shows a 
modulation amplitude consistent with zero in models lb and 2a. 

6.2. All bins analysis: dark matter or other physics? 

When the data from all three energy bins are analysed simultaneously, we find that only 
the DM models la and lb receive weak support from the data over no modulation. In 
contrast, models 2a and 2b are moderately disfavoured with respect to the no modulation 
scenario. DM models are thus strongly favoured over "other physics" models, with odds 
in favour of the DM models ranging from 185 : 1 to 560 : 1, see table 6. This is a 
consequence of the predictiveness of models la and lb; that is, Occam's razor is at work, 
penalising the "other physics" models for their excessive free parameters unsupported 
by the data. Table 6 summarises the odds for each model against the no modulation 
model 0. Odds between any pair of models can be obtain simply by multiplication of 
the odds given there. 

The inferred values for 5**^ are closely in line with the single bin analysis (see 
figure 5), with the exception of bin 3 in model 2b, which now prefers = 0.03 ± 0.04 
in the combined fit, instead of Sf^ = 0.22 ± 0.04 when the bin is analysed on its own. 
This result comes about because the phase tmax and the period T are constrained to be 
the same in all three bins in the combined fit. The fit, in turn, is driven mainly by the 
data in bins 1 and 2. Indeed, the common phase preferred by the combined data is now 
^max = 104 ± 10 days for model 2a and 107 ± 12 days for model 2b, while the preferred 
period is T = 344 ± 22 days in model 2b (see also figure 6) . Figure 7 also shows that 
the all-bin fit for tmax and T is dominated by bin 2, which singles out a value tmax ~ 100 
and T ~ 345 from the vast degenerate region obtained from bin 1 alone. On the other 
hand, the multi-modal posterior preferring T ~ 100 days from bin 3 is completely cut 
away in the combined fit. This shows that bin 3 on its own prefers value for phase and 
period that are incompatible with bins 1 and 2. 

It is important, like in any good Bayesian analysis, to assess the robustness of the 
results with respect to reasonable changes in our prior choices. An interesting question 
to ask in this regard is, have the Bayes factors for models 2a and 2b been artificially 
suppressed because of our choice of priors? To answer this question, we appeal to 
the Savage-Dickey density ratio defined in equation (2.6). Consider for concreteness 




Figure 7. 2D marginal posterior pdfs in the {imax, T'j-plane for model 2b from the 
single-bin analyses — bin 1 (top left), bin 2 (top right), and bin 3 (bottom left) — and 
for the combined fit (bottom right). The black solid lines are 90%-contours. 



model 2a. If we were to reduce the prior ranges for all three fractional modulation 
amplitudes to 5*^ = —t- 0.5 from the default choice of 5*^^ = — ?■ 1, then it follows from 
the SDDR formula that the Bayes factor In B in favour of model 2a would increase by 
approximately In 2'^ ^2.1 units, bringing it to ~ —1.06 relative to model 0. In such a 
case, the model comparison between model 2a and the consistent DM model lb would 
produce a weak evidence in favour of the latter scenario, while a comparison between 
models 2a and la would still favour moderately the latter. Indeed, even with the prior 
ranges reduced further to 5*^ = — )■ 0.2, model 2a would only barely overtake model lb, 
and would still lag marginally behind model la as the most preferred model by the data. 
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Table 6. Odds between modulation models A4i versus the no modulation model A4q 
for the CoGeNT modulation data. 



Model Mi 



Mi -.Mo 
Bin 1 Bin 2 Bin 3 All 3 bins 



la 
lb 



2:1 4:1 1:1 8:1 

2:1 2:1 1:1 5:1 



2a 
2b 



1:7 1:1 1 : 16 1 : 37 
1:9 1:3 1:6 1 : 70 



Therefore, we conclude that the general statement that DM-like models are preferred 
over "other physics" models is robust from a Bayesian point of view. 

6.3. Comparison with classical p-values 

In frequentist statistics, classical hypothesis testing seeks to rule out the null hypothesis 
Hq by quantifying the probability of observing data as extreme or more extreme than 
what has been obtained. To this end, a test statistic t can be defined, in such a way that 
extreme values of t are increasingly improbable under Hq. For concreteness, we assume 
that, under the null, larger values of t have monotonically decreasing probabilities. Then 
the tail probability, or the p- value p, can be computed via 



so that small values of p denote that the observed data are very improbable under the 
null. We stress once more that p-values are not probabilities for hypotheses — they are 
probabilities of obtaining more extreme data than observed assuming the null hypothesis 
is correct. In order to obtain the probability for a hypothesis (which is arguably the 
scientific question we are interested in), one needs to take a Bayesian approach, as we 
do in this work. Also, we caution that the mapping of the test statistic t onto a p-value 
requires in general a Monte Carlo simulation; analytic solutions exist only in special 
cases, and apply only if certain regularity conditions hold, as we discuss below. 

A popular choice for the test statistic t in the context of nested models is the 
effective chi-square (or more precisely, the profile likelihood ratio), || 



unconditional maximum likelihood in the whole parameter space of the more complex 
model. Let us identify the simpler model (with i? = i?*) as the null hypothesis that 

II Strictly speaking, the quantity we compute for model lb is the posterior ratio, since the prior pdfs 
in that case are not uniform. 




(6.1) 
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Table 7. '^xle values, as defined in equation (6.2), for the various modulation 
scenarios relative to model 0. Where analytically calculable, we quote also the 
corresponding classical p- values of the null hypothesis (model 0), obtained via 
Chernoff's theorem, Eq. 6.3. We also give the number of extra free parameters in 
the alternative hypothesis relative to the null. 







Axgjj relative to model 


Model 


Bin 1 


Bin 2 


Bin 3 


All 3 bins 


la 


2.04 


4.23 




6.26 




p = 0.08 p = 0.02 




p = 0.02 




(^ = 1) 


(^ = 1) 




(^ = 2) 


lb 


1.94 


1.88 


0.020 


3.84 




p = 0.0< 


B p = 0.09 


p = OA 


p = 0.1 




{u=l) 


{u = l) 


{u = l) 


{u = 3) 


2a 


3.61 


8.36 


0.025 


10.63 


2b 


3.70 


8.87 


10.88 


10.86 



we seek to rule out (e.g., no modulation in the CoGeNT data). If the likehhood in 
the N additional parameters of the more complex model is Gaussian and unbounded, 
then Wilks' theorem [69] applies, meaning that the test statistic Ax^g is asymptotically 
distributed as a with N degrees of freedom. We note however that one of the 
conditions that validate the application of Wilks' theorem is that the likelihood must 
be unbounded, i.e., the additional parameters of the more complex model cannot sit on 
the boundary of its parameter space. This is precisely the situation we encounter when 
treating the modulation amplitudes 5**^, where 5*^^ = under the null hypothesis. Hence, 
Wilks' theorem cannot be applied to obtain the p-valuc from the Axgg, a procedure 
that is often followed outside its realm of vahdity (see [70] for a discussion aimed at 
astronomers) . 

For the case where (i) the N additional parameters are bounded, with the 
null hypothesis sitting on the boundary, (ii) the likelihood is Gaussian, and (iii) all 
parameters arc identifiable under the null hypothesis, we can use instead Chernoff's 
theorem [71,72] to evaluate the p-value of the null. Chernoff's theorem says that the 
distribution of the test statistic Axlg under the null is asymptotically a weighted sum 
of random variables xl following chi-squared distributions with i degrees of freedom: 

P = E2-^f^)p(x^>AXe\)- (6-3) 

Equation (6.3) can be used to compute the p-value of the null hypothesis of no 
modulation when the more complex hypothesis is identified with models la or lb, i.e., the 
DM hypotheses. However, it cannot be applied to the "other physics" models 2a and 2b, 
for these models contain parameters (the phase and the period) that are undefined and 
unidentifiable under the null, i.e., when Sl^ = in models 2a and 2b, the parameters 
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tmax and T are meaningless. 

Applying equation (6.3), we obtain a p- value of 0.02 under the null for the three 
energy bins combined, when the alternative hypothesis is model la. This formally 
corresponds to a 2.3a detection (using a Gaussian distribution to convert p-values into 
the number of sigmas). If instead we take model lb as the alternative, the p- value is 0.1, 
equivalent to a l.Gcr detection. Monte Carlo simulations would be required to determine 
the distribution of the test statistic when the alternative model is either 2a or 2b. 



6.4- Impact of an anisotropic velocity distribution 

Although DM models appear to be preferred over modulation models due to other 
physical processes, there remains some conflict between the modulation amplitude 
preferred by the CoGeNT time-dependent data in energy bin 2 and the amphtude 
predicted by the experiment's total event rate. Therefore, we wish to investigate 
whether there is a WIMP scenario that could reconcile these findings. Modifications 
to the particle physics interaction have been widely discussed in the context of, e.g., 
inelastic or isospin violating DM models [15-21,30]. Here, we consider the possibility 
of an anisotropic DM velocity distribution in the Galactic halo. Anisotropic velocity 
distributions have been investigated in [13] as a means to reconcile the CoGeNT and 
the DAMA data, and more recently in [32] to reconcile both the CoGeNT total rate 
with the CRESST excess and the CoGeNT modulation with its total rate. 

Consider a class of ellipsoidal (or equivalently, triaxial) halo models whose 
gravitational potential is given by ^{x,y,z) — {vc/2)\n{x'^ + y'^/p^ + z^/q^), where 
X, y, z are the (Cartesian) coordinates, and < p^ < 1 parameterise the degree 
of triaxiality [73]. Assuming the Earth sits on either the long axis (x-axis) or the 
intermediate (y-axis), the corresponding DM velocity distribution is in the form of a 
triaxial Gaussian, i.e.. 



/(^'W) = 7TT-Ti7^7— — exp 



2al 2al 2al 



(6.4) 



in the Earth's rest frame, where is the Earth's speed in the Galactic frame given in 
equation (4.2), and ax-, with X = R,(f),z, is the velocity dispersion in the X-direction 
in a cylindrical coordinate system. The exact expressions for ctr ^^^ can be found in 
equations (4.25) and (4.26) of [73]: specifically, they depend on vq and on the anisotropy 
parameter 7. In the spherical limit, equation (6.4) reduces to the well-known Gaussian 
velocity distribution of singular isothermal halos.^ 

One possibility is a radial anisotropy, ctr > cr^ — crtf>, or equivalently, 7 < 0. 
However, this kind of anisotropy does not give the desired effect of enhancing the 
fractional modulation amplitude, because the WIMP phase space is partially depleted 
in the direction of the Sun's movement around the Galactic centre, thereby reducing 

^ It is a reasonable approximation to consider an isothermal halo and its ellipsoidal extensions, because 
direct detection is sensitive only to the local properties of the velocity distribution and cannot presently 
distinguish between different DM density profiles [45, 50]. 
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Figure 8. 3D marginal posteriors for {todm, S'^'^o/^}, where the S^^ direction is 
indicated by the colour coding, for a specific realisation of an anisotropic DM velocity 
distribution (see main text for details). The left panel shows S^, and the right S"^^. 
The arrow in each plot points to a circled region in which the predicted values of 
and Sf^ find good agreement with the CoGeNT modulation data. 

the fractional difference between the maximum and the minimum event rates relative 
to the mean [74]. In contrast, a tangential anisotropy 7 > may potentially enhance 
the fractional modulation amplitude [52, 73, 74]. 

We perform a MCMC for one specific realisation of a tangentially anisotropic 
triaxial halo (anisotropy parameters p = 0.9, q = 0.8 and 7 = 16, and astrophysical 
parameters fixed at their mean values, see table 4) to infer the joint posterior pdf of 
the DM mass and cross-section — as well as the nuisance parameters — preferred by the 
CoGeNT total event rate. The DM mass is sampled using a uniform prior in the range 
1 — 7- 20 GeV, while a uniform prior on log(o"^Ycm^) in the range —42 — )■ —38 is assumed 
for the spin-independent interaction. Along the lines of model lb, we derive posterior 
pdfs for the fractional modulation amplitudes S"^ via equation (5.3). 

Figure 8 shows the resulting 3D marginal posteriors for {rriDM, c"n^ ^m} ^'^^ ^^^t 
two energy bins, AEi and AE2. Clearly, the predicted fractional modulation amplitudes 
in energy bins 1 and 2 have increased has a result of the anisotropic velocity distribution. 
In bin 1, the modulation amplitude preferred by the total rate is = 0.14 ± 0.03 
(compared with S"^ = 0.098 ± 0.021 for the isotropic DM model lb), slightly larger 
than the ~ 10 — )■ 11% modulation favoured by the CoGeNT time-dependent data (see 
figure 5). For energy bin 2, we find S*^ = 0.05 ± 0.02, which is an improvement on 

= 0.026 ± 0.011 predicted by model lb, but still significantly lower than the 20% 
required to explain the CoGeNT time modulation. In bin 3, = (0.12 ± 0.77) x 10^^ 
is consistent with no modulation and also with the expectations of models la and lb. 

Interestingly, although the "best-fit" S^^ and especially values are not perfect 
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matches to the fractional modulation amplitudes required to fit the CoGeNT modulation 
data, we see from the 3D marginal posterior for {tbdM) o"n^ ^m} figure 8 that it 
is possible to find points in the preferred region in the DM parameter space that 
could improve further the agreement between the CoGeNT total and modulated rates. 
Take as an example mDM — 10 GeV and ~ 10"'*° cm^. These numbers lead to 
~ 0.10 0.13 and 5^ ~ 0.12 0.15 (see the circled region in figure 8). Prom 
this, we conclude that an anisotropic DM velocity distribution appears to be a promising 
direction towards reducing the tension between the CoGeNT total and modulated rates. 
We defer a more complete exploration of this avenue to a future work. 

7. Conclusions 

The annual modulation signal claimed by the CoGeNT collaboration is an intriguing 
puzzle for the DM direct detection community. In this paper we have investigated 
the origin of the time-dependent signal detected by CoGeNT, and performed a 
sophisticated model comparison analysis to identify the best physical explanation. Our 
phenomenological approach features three nested scenarios: no modulation (amplitude 
of modulation set to zero), modulation due to DM (phase, modulation amplitude and 
period predicted by WIMP models), and modulation due to other physics (amplitude, 
phase and period of modulation as free parameters). We have used the Baycsian model 
comparison framework to assess the strength of evidence in favour of these scenarios. 
In addition, we have recomputed classical p-values against the null hypothesis of no 
modulation where analytically possible, since they had been misestimated in previous 
works. 

Regarding the evidence for an annual modulation with a maximum rate occurring 
close to June 2 (or t 

max ~ 152 days), we find that there is weak evidence for a modulation 
in both the energy ranges 0.5 0.9 keVee and 0.9 — > 3.0 keVee. This conclusion is borne 
out both from a Bayesian model comparison point of view and from a classical hypothesis 
testing perspective. The energy range 3.0 — )■ 4.5 kcVcc shows no sign of modulation with 
a period of one year. Modulation models due to "other physics" compare unfavourably 
with the no modulation case, paying the price for their excessive complexity. 

The question of whether this signal is due to DM or to other physical processes 
is more difficult to assess. When all energy bins (0.5 — > 4.5 keVee) are considered, 
the DM models are strongly favoured over the phenomenological models in which the 
signal is due to other physics, with odds of several hundreds to one, depending on the 
details. This conclusion is fairly robust with respect to changes in the priors assigned 
to the parameters of the "other physics" models, in the sense that reasonable changes 
to the prior ranges could at most reduce the preference but not overturn the ranking. 
Interestingly, between the "naive" DM model la and the "consistent" model lb which 
enforces self-consistency between the modulation signal and the CoGeNT total rate, 
there is no compelling preference for the latter, despite the fact that compatibility with 
the total rate dictates that the modulation amplitude be significantly smaller than the 



Evidence for dark matter modulation in CoGeNT? 



25 



observed amplitude in the energy range 0.9 3.0 keVee. Nonetheless, we have shown 
that this tension between the total and the modulated rates could be partially remedied 
by using a triaxial DM velocity distribution with a tangential anisotropy. 
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